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ABSTRACT 


Under an age replacement policy a system is replaced at a fixed age <p or at failure 
whichever comes first. If the cost of replacing the system before failure is less than the 
cost of replacing it at failure, this type of maintenance policy can lead to considerable 
savings. An often used criterion for finding an "optimal" replacement age (p, is to mini¬ 
mize the long run expected cost per unit time of a policy with replacement age cp. This 
cost function clearly depends on the underlying distribution of the system lifetimes. 
When this distribution is unknown, the cost function and hence (p* need to be estimated. 

In this thesis, we study the large and small sample properties of a procedure which 
estimates cp*. In particular, we study sequential maximum likelihood estimators of (p* 
which are updated at each replaccm''nt based on the replacement history of the system 
so far. In this sequential procedure each system is subject to the age replacement policy 
with estimated (p* based on all the data gathered so far. This type of procedure should 
control the actual cost per unit time while gathering data needed to estimate (p*. 

This thesis contains a detailed description of the sequential estimation procedure 
when the underhing system life times have a Weibull distribution and a Gamma dis¬ 
tribution. Monte-Carlo methods are then used to study the behavior of the estimated 
optimal age replacement times and more importantly the actual costs per unit time for 
different sample sizes, costs and choices of the underhing Weibull and Gamma distrib¬ 
utions. 
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I. INTRODUCTION 


A. BACKGROUND 

Optimal maintenance policies are designed to reduce the number of system failures 
and minimize the cost of repair by scheduling planned replacements. By far, most of the 
research in this area has been from the modeling stand-point. Even in the most basic 
scenario, where the underlying system lifetimes are assumed to be independent and 
identically distributed, the problem of updating the maintenance policy using the past 
maintenance history has not been adequately solved. In this thesis, Monte-Carlo 
methods are used to study a particular parametric procedure which updates estimates 
of an optimal age replacement policy after each replacement. 

B. MAINTENANCE AND MAINTENANCE POLICIES 

Maintenance is a combination of actions carried out to retain a unit in, or restore 
it to. an acceptable condition. There are two forms of maintenance. These are prevemive 
and corrcciivc maintenance. Corrective maintenance is carried cut to restore (including 
minor adjustment and repair) a unit which has ceased to meet an acceptable condition. 
Preventive maintenance is interned to reduce the likelihood of a unit not meeting an 
acceptable condition [Ref 1: pp. 1-8]. 

Maintenance policies reduce the number of system failures and maintenance costs 
by adopting a schedule of planned maintenance. For instance, when the failure of a unit 
during actual operation is dangerous or costly (i.e., the cost C, of unscheduled mainte¬ 
nance due to system failure is more than the cost Cj of scheduled maintenance before 
system failure) and the unit is characterized by a failure rate that increases with age. it 
may be wise to maintain it before it has aged too greatly [Ref 2: p. 46]. On the other 
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hand, performing maintenance actions too frequently can also be costly. Thus, choosing 
the preventive maintenance policy that schedules maintenance to best control costs will 
be influenced by the relative costs of failure and preventive maintenance and the 
stochastic properties of the lifetime of the unit [Ref 3: pp. 19-24]. We are interested 
in determining the sequence of times at which preventive maintenance should occur. In 
particular, we are actually interested in determining an optimal maintenance policies 
which minimizes long run expected maintenance cost per unit time. 

Preventive maintenance is the total of all service functions aimed at maintaining and 
improving reliability performance characteristics and concerns itself with such activities 
as the replacement and repair of systems, inspections, testing and checking of working 
parts during their operation. In this thesis, we will only consider maintenance actions 
which involve replacement of a system. It will be assumed that the replacement action 
returns the equipment to the as ne>\ condition, thus providing the same services as the 
equipment which has just been replaced. By making this assumption, we are implying 
that the distribution of time to failure of the new s}stcm is the same as that of the system 
which was replaced. In addition, we assume that system lifetimer are independent. The 
unscheduled and scheduled replacement costs C, and Cj remain constant where C, > Cj. 

.An important replacement policy is the policy based on age {age rcplaccmem). Such 
a polic) is in force if a unit is always replaced at the time of failure or cp units of time 
after its installation, wliichever comes first. Under a block replacement policy the unit 
is replaced at times k(p (k = 1,2,...), and at failure. This replacement policy derives its 
name from the commonly employed practice of replacing a block or group of units in a 
system at prescribed times k<p (k = 1,2, ...) independent of the failure h’Story of the 
system [Ref 2: p. 46]. 
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'' Age replacement is administratively more difficult to implement, since the age of 
the unit must be recorded. But block replacement, although simpler to administer since 
the age of the system need not be recorded, leads to more frequent replacement of rela¬ 
tively new items " (Ref. 4: p. 158]. In this thesis we only consider age replacement 
policies. 

An optimal age replacement policy is the age replacement policy which yields the 
sm.allcst long run expected replacement costs per unit time. To find an optimal age re¬ 
placement policy we require explicit knowledge about the system's lifetime distribution. 
When we don't know the lifetime distribution explicitly, the systems life distribution and 
the optimal age replacement policy needs to be estimated. Estimation based on a fixed 
sample of independent and identically distributed {i.i.d.) system life times has been ex¬ 
amined in detail [Refs. 5,6,7.8]. In order to collect such i.i.d. data, experimental systems 
must be left in service until failure. When observation of complete system lifetimes is 
not available (because it either takes too much time or is too costly), the most cost ef¬ 
fective approach is to start with an initial estimate of the optimal replacement age and 
then update this estimate after each system replacement. After each replacement, the 
next s}stem is subject to a replacement polic\ that is close to the best estimated polic\ 
so far (Ref 9: pp. 2-3], This procedure is described in detail in Chapter 11. We will use 
simulation to stuJ\ this sequential estimation procedure when the undcrljing life dis¬ 
tribution comes from a parametric family. In particular, in Chapter 111 we study this 
procedure with an underlying Weibull distributions and in Chapter IV we consider 
Gamma life distributions. Conclusions and recommendations are given in Chapter V. 
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II. THE SEQUENTIAL ESTIMATION PROCEDURE 


A. THE OPTIMAL REPLACEMENT AGE 

Consider a system or component which is replaced upon failure. When the cost of 
a replacement that is planned in advance is less than the cost of an unplanned replace¬ 
ment, a simple age replacement policy can lead to considerable savings. In this type of 
policy, an age cp is specified; items that are still functioning at that age are replaced 
(these are planned replacements); items which fail and are thus replaced prior to are 
the unplanned replacements. 

Let A'j,... be a sequence of independent and identically distributed {Lid.) positive 
random variables with common distribution function Fg belonging to a family of dis¬ 
tributions parameterized by 0. The sequence Aj,... represents the sequence of system 
lifetimes that would be observable if the systems were replaced at failure. Let C„ C 2 . ( 
C, > C 2 ) be the respecti\e costs of planned replacement (before system failure) and un¬ 
planned replacement (at sjstem failure). The observed durations between replacement, 

min(A'. c?) i = 1, 2.form a renewal processes [Ref 2: p. 87]. Therefore, the long 

run expected cost per unit time with age replacement at cp is 

C, X Foiw) -f C2 X Fff((p) 

CUp) = --- 

El min (A), (p) ] 

^ C,xF,i<p)+C2xF,{(p) ( 2 . 1 ) 

dx 

0 

where Fe{x) - 1 — Fg{x) is the survival function. In equation (2.1) the numerator is the 
expected cost of one replacement under the age replacement policy, and the denominator 
is the expected time between replacement. 
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Under reasonable conditions, there exists a unique and finite time cp* <oo where 
C{(p) attains a global minimum (Ref. 10:, pp. 161-168]. For example, a sufficient con¬ 
dition for the existence of (p* is that has a failure rate /(r) that strictly increases to 
infinity. 

B. THE SEQUENTIAL ESTIMATION PROCEDURE 

Finding the optimal planned replacement age cp* requires knowledge of the under¬ 
lying distribution function Fg. If the underlying lifetime distribution Fg is completely 
specified, then finding <p* can be done either explicitly or numerically using equation 
(2.1). However, if the parameter 6 is unknown, we need to estimate (p*. If full lifetimes 
A'l, Aj, .... X„ are available, then C{(p) and thus (p* can be estimated by replacing 0 in 
equation (2.1) with the usual parametric estimators of 6. This approach has the obvious 
disadvantage that it requires several s} stems to operate until failure to be observed be¬ 
fore a replacement policy is implemented. In addition, subsequent data is not used to 
update the policy. A more practical and potentially more cost effectit e approach is to 
update estimate (p* after each replacement and implement the updated age replacement 
policy (through estimates of cp*) after each replacement. 

Let (Pi, (p,, ... be a sequence of estimators of (p where depends on the first n re¬ 
placement ages. A procedure to compute the estimators is developed as follows: 

1. At the «th replacement observe the system lifetime X„ or whichever comes first. 
Let Z„ = min(A'„, and d„ = /(A'„ < where 1(A) is an indicator function of 
the set A. In other words, if the unit is replaced before failure, then the replace¬ 
ment time Z„ = <p„., and = 1, otherwise Z„ = A' and <5„ = 0. 

2. The data available to estimate are the pairs (Z„S,) i= 1, 2, ..., n (Throughout, 

without loss of generality we take special case as Z, = A',, = I). The maximum 

likelihood estimator (MLEj 6„ of 0 is computed from this right censored data. 

3. is then found to minimise equation (2.1) with 0„ taking the place of the unknown 
parameter 0. 
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The procedure is then repeated. The technical problem of using this data to estimate cp 
at each stage is that the pairs (Z„ 6,) i= 1, 2,.... n are clearly not i.i.d. Thus, the usual 
methods for studying the properties of the sequence of estimators {(p„} from the right 
censored data are not immediately applicable. Thus, we will use simulation to study 
both the large and small sample properties of this sequential estimation procedure. 

The replacement cost for the /th system is C, if X, < ^,.i, otherwise the replacement 
cost is Cj. With this sequential estimation procedure the actual total replacement cost 
for the first n systems that are observed is given in equation (2.2) 

n 

+ ^2 X ( 1 - ) }. (2-2) 

1=1 

and equation (2.3) is the total operating time for the n systems. 


n n 

(=I •■=! 

When studying the properties of this sequential estimation procedure it is important to 
sec if and how fast the actual cost per unit time -y- converges to optimal cost C{(p*} 
(Ref llj. 

In this thesis we simulate the sequential estimation procedure for two parametric 
families of distributions considered: In Chapter III, F belongs to a Wcibull family of 
distributions with the shape parameter a > 1 and unknown scale parameter (/), and in 
Chapter IV, F belongs to a Gai,ima famih of distributions with shape parameter p > 1 
and unknown scale parameter (0). Both the Weibull and Gamma distributions were 
chosen for the simulation because they arc commonl) used to model ssstem lifetimes. 
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they have increasing failure rates when their shape parameters > 1, and because esti¬ 
mation of the unknown parameters and minimization of the estimated cost function are 
numerically tractable. 
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III. SIMULATION SETTING FOR THE WEIBULL DISTRIBUTION 


A. UNDERLYING LIFE DISTRIBUTION 

This chapter is concerned with the estimation of the optimal replacement time when 
it is known that the underlying lifetime distribution is a member of the tw'o parameter 
Weibull family with shape parameter a and scale parameter ). > 0, w’here the density is 
given by 

f(i) = a j. e~ for t>0. (3.1) 

The Weibull distribution has failure rate 

/(/) = a?. , r>0. (3.2) 

When a > 1.0, the failure rate in equation (3.2) is strictly increasing to infinity. Thus, 
for Weibull distributions, with a > 1.0, a unique and finite optimal replacement age cp* 
exists. To compare the results with [Ref 9; p. 12], the same ten different Weibull dis¬ 
tributions used in the simulation have a values 1.1, 1.2, ..., 1.9, 2.0. This selection of a 
values gives us a range of distributions which become more like the exponential distrib¬ 
ution as a decreases from 2.0 to 1.1. The Weibull distribution with shape parameter 
a = 2.0 is called the Rasleigh distribution and has a linearly increasing failure rate. To 
make fair comparisons between Weibull distributions, the scale parameter / was chosen 
so that the expected system lifetime E(.V,) = 2.0. In our figures we have selected only 
certain values (a = 1.2, 1.4,1.6,1.8, 2.0) for the scale parameter a. This was done to 
make the figures more readable. See Figures 1 and 2, for plots of the Weibull densities 
and corresponding failure rates. 
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TfEIBULL DENSITY FUNCTION 
f(t) = 



Figure 1. The Weibull Density Function f(t) >uth £(A’) = 2.0 
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WEIBULL DENSITY FAILURE RATE FUNCTION 


\[t) = 



Figure 2. The Failure Rate of The Weibull Distribution with = 2.0 


B. OPTIMAL REPLACEMENT TIME 

When the distribution of the system lifetime is Weibull, the reliability (survival) 
function defined by 


i>0, ;. >0 , ( 3 . 3 ) 

1 ^ 0 . 

The long run expected cost per unit time under a simple age replacement policy with 
scheduled replacement at age (p, by using equation (2.1) is given by 


F{t) = P{X>i) 
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Q,» = 


C,(l - + 


r 

‘'n 




(3.4) 


i/x 


See Figure 3, for plot of the long run expected average cost function when the underly¬ 
ing life distribution is Weibull with shape parameter a varying from 1.1 to 2.0. For each 
curve on Figure 3 the optimal replacement time cp* can be located on the x-axis at the 
minimum point. 



Figure 3. The Long Run Expected Average Cost Curves with £(A’) = 2.0 
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C. SEQUENTIAL ESTIMATION PROCEDURE OF THE WEIBULL 
DISTRIBUTION 

We will use the sequential estimation procedure which we discussed in general in 
Chapter II, with the known shape parameter a > 1.0 and the unknown scale parameter 


Let be the MLE of X after n replacements. The corresponding estimator of 
(p* which minimizes C; J^cp) needs to be found numerically. The following result, sim¬ 
plifies this procedure considerably. 

Lemma 7. Let (p* minimize then (p* = Xx* where x* minimizes Ci,,(jr). 

PROOF, By using equation (2.1), the cost function at age t with scale parameter / be¬ 
comes. 


Q ,=,(0 = 


f' 

/{x) dx -I- Cj Ax )' 

_ 

r/ Poo 

y(v) dv du 

Vfj 


If wc substitute the Weibull density into the equation (3.5), the numerator can then be 
shown to be 


C] /’ax* dx + C 2 X.^ax^~^e dx. 


with the change of variable y = /a-, wc obtain 
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(3.6) 




If we factor the constant term ( — )•’"* from (3.6), we get 


In expression (3.7) the constant terms cancel, and the integrands are just the density of 
the Weibull distribution with / = 1. Thus (3.7) is 

C,f,=,(;.0 + Qf,=,(//). (3.8) 

With the change of variables y = /v the integrand in the denominator of equation (3.5) 
becomes, 


ay“ dy = /•;=,(/r). 

?.U 

From (3.8) and (3.9), the cost function may be written as 


Q ..(0 


Q ^.)= i (^.0 + Q ^;.= i (^-0 

P' 

F,^^{/.u)du 

•^0 


(3.9) 


Set y = /u, then 
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(3.10) 





= ^.C.= 





Therefore, if x* minimizes, and <p* minimizes Cj,,(r) then (p* = ?.x*, ■ 

To estimate (unknown scale parameter) from the right censored data; recall that, 

A 

the available data to estimate are the pairs of (Z„ ^,) i= 1, 2, ..., n where Z, = and 
= 0. By using the maximum likelihood estimation (MLE) procedure we can estimate 

A 

/, from the first observation A',. The MLE of /. can be found by differentiating the 
likelihood (or the loglikelihood) with respect to /, setting the results equal to zero, and 
solving for /. 

L^;.,a|Z) = 


/(/, a\Z) = In a + a In >. + (a — 1) In Xj — /“x*, 


cl{).. a\Z) 
ck 



1 

A' 


a 
1J 


c/(/,a|Z) ^ A I 

--- = 0 => /, = —. 

cV. ‘ -^1 

By Lemma 1, we can estimate (p* such that = >.,x* where minimizes G (cp) 
and X* minimizes C, ,(x)., 

After the second observation, we set Zj = min(A 2 , ^*) and 

(• 1 A2<^>r 

(. 0 otherwise. 
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After n observations, data will be (Z,, (5,), (Zj, (5j). (Z„, where Z, = min(A'’, , 

<p*,), Z, is the /til lifetime, and (p*^ is the best estimation of our optimal policy so far. 


In general, if we repeat the sequential estimation procedure n times, we obtain 
LiX, a|Z) = /( 2 ,) [ 7 ( 22 )]'^ ... 

l{X.,a\Z) = «(1+(52 + ... + ^„) In a — (z“ + 22 + ... +z^), 




By Lenuna 1, we can find the MLE of cp* explicitly from the MLE /„ of ). based 

on the first n replacements using 

(p„ = A„xx . (3.12) 
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D. SIMULATION RESULTS FOR THE WEIBULL DISTRIBUFION 


1. Finding of The Scale Parameter ), 

As mentioned before, the ten different Weibull distributions used in the simu¬ 
lation have a values 1.1, 1.2, ..., 1.9, 2.0 and the scale parameter X was chosen so that 
the expected system lifetimes £(A') = 2.0. The mean of the Weibull distribution is 


£CA'] = 


r(T) 


(3.13) 


In equation (3.13), the scale parameter X value can be obtained easily since the scale 
parameter a and £CA'] are known. The APL program Weibull in Appendix A calculates 
the scale parameter X values for given a and The results are given in Table 1 for 

the shape parameters and corresponding scale parameters. 


Table 1. SCALE PARAMETER / VALUES FOR £(A:,) = 2.0 


Shape pa¬ 
rameter a 

a 2.0 

a= 1.9 

a= 1.8 

a = 1.7 

a = 1.6 

Scale pa¬ 
rameter / 

0.44311346 

0.44368166 

0.44464337 

0.44612225 

0.44828714 

Shape pa¬ 
rameter a 

a = 1.5 

a = 1.4 

a= 1.3 

i 

a= 1.2 

m 

Scale pa¬ 
rameter X 

0.45137265 

0.45571117 

0.46178836 

0.47032793 

0.48245624 
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2. Finding of The Actual Optimal Age Replacement Time and Corresponding Actual 


Cost 

The actual age replacement time (p* can be located in Figure 3 on the x-axis 
at the minimum point of the cost function Q,.(<p) (3.4). In order to find the unique 
minimum point cp*, it hard to solve equation 3.4 analytically. Therefore, the actual age 
replacement time (p* found by simulating the cost function (3.4) for different (p where 
the simulated cost function is given by 


C, X 




+ - ] 


r 


^min(X,-, cp) 


(3.14) 


where X, (i = 1, 2. ..., n) are the simulated i.i.d. Weibull random variables with the spe¬ 
cific a and /. cp* is the minimum of Q,,(^). Since the optimal replacement time cp* 
comes from the simulation result, it varies slightly with the number of pseudo random 
variables used and the seed numbers used to generate them, thus we choose 
« = 15 X 10-'. 

3. Finding of The Average Age Replacement Times and The Average Costs 

In order to find the average of the estimated optimal age replacement times and 
the average actual costs based on lifetimes with a Weibull distribution, wc wrote the 
Fortran simulation program Averweib in Appendix C. From Lemma 1, wc know that 
X* minimizes C^^,_,(/) and (p* minimizes Q,(/) when ip* = }.x*. Thus, x* can be found 
dividing the actual optimal age replacement times by the scale parameter z, where the 
actual age replacement time is taken from the simulation program Sim and the / values 





are taken from the Table 1 for the specific value of the shape parameter a. The un¬ 
planned and planned replacement costs C, and Q (to calculate the cost function); the 
scale parameter / (to calculate the x* value); the actual age replacement time (from 
program Sim to calculate the x* value and to find MSE for each estimated age replace¬ 
ment time); the actual cost value (from program Sim to find MSE value for each esti¬ 
mated cost); and the shape parameter a must be given by the user in the initialization 
part of the simulation program Aweru'eib. In much of the simulation, we considered the 
sample sizes A' = 10, 50, 250 and 1000. Other sample sizes are also considered, but in 
much less detail. Each repetition of the simulation is based on generating A' system 
lifetimes. These system lifetimes (A') are used one at a time for the sequential estimation 
procedure. At the first observation, Z, = A\, and (5, = 1. The unknown scale parameter 

A 

/ values are estimated by using the equation (3.11). After finding /, the estimated age 
replacement time values are calculated by using the equation (3.12) in the simulation. 
By using A’ generated system lifetimes, we deterntine A' estimated age replacement times 
and A' estimated costs in each simulation. We repeat this simulation 1000 times 
(NREP= 1000) and then we find the average value for both age replacement times and 
replacement costs. 

Let C^,v, j = 1, 2, ..., 1000 be the actual cost per unit time for the first A' re¬ 
placements of the jth repetition of a simulation (where each Qv is computed using (2.2) 
and (2.3)). The average actual cost per unit time over the 1000 repetitions is 


C,v = 



(3.15) 
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Let <p*,j = 1,2. 1000 be the estimated optimal replacement time for the first 

N replacements of the yth repetition of a simulation. The average age replacement time 
over the 1000 repetitions is 


1000 



y=i 


For each simulation, we also calculated 


1000 / 

Z (C;v-C((p )) 

^-iooo-• 


y=i 


(3.16) 


(3.17) 


.\1SL.\GE 


1000 , * *k2 


E 

;=i 


1000 


(3.18) 


where .MSLCOS is the a\cragc squared difference of the actual replacement cost per unit 
time from the estimated minimum long run expected replacement cost per unit time and 
MSLAGL is the average squared difference of the actual optimal age replacement time 
from the estimated optimal age replacement time. These MSB values are calculated in 
the simulation in order to see if the sequential estimation procedure is converging the 
actual cost and the actual optimal age replacement time. If we get the MSB values close 
to the zero, then we can say that this procedure is working well. 

Tables 2, 3, 4 and 5 sununarize the simulation results of the optimal age re¬ 
placement times and the minimum long run expected replacement costs per unit time for 
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difTercnt values of a with N= 1000 for fixed costs C, = 2.0, C, = 5.0, C, = 8.0, C, = 10.0 
while holding the Cj fixed at 1.0. Tables 6, 7 and 8 summarize the simulation results 
using the three sample sizes of 10, 50 and 250 as small, moderate and large sample sizes 
for different values of a with fixed costs C, = 2.0 and Q = 1.0. Included in Tables 
2~8, is the probability that a system will fail before the time (p under the optimal age 
replacement policy. 


P{Xi<(p*) = 1 >*. (3.19) 

We have also plotted the results of the average age replacement times (Fig. 4) 
and the average costs (Fig. 5) for a= 1.2, 1.4, 1.6, 1.8, 2.0 when fixed costs C, = 5.0 and 
Cj = 1.0. From these plots and the results from Tables 2~8, we observe that in general 
when a decreases from 2.0 (i.e., the underhing life distribution is becomes more expo¬ 
nential) the optimal age replacement time cp*, the long run expected optimal replace¬ 
ment cost and .VISE values increase. For values of a. close to 1.0. very little is gained by 
replacing the system before failure at the higher cost C,. When the ratio of the un¬ 
planned replacement cost C, to the planned replacement cost Cj increases the optimal 
replacement time for a specific a decreases. The larger values of <p*, insure that a small 
percentage of replacement will be made before failure, which is what we desire if the 
system's life distribution is close to exponential. 

By looking at the tables for different sample sizes, we also observe that the av¬ 
erage cost per unit time up to the A'th replacement decreases with A', the number of re¬ 
placement. This result is promising because the ultimate goal of the sequential 
estimation procedure is to decrease costs while sampling. Even though as A’ -* oo, the 
long run average co.si per unit time will approach the optimal replacement cost C{(p*) 
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with probability 1.0 [Ref. 11], there is no guarantee that the average replacement cost 
will decrease for the first observations. 
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Table 2. 


ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
WEIBULL MODEL WITH C, = 2.0, C, = 1.0, E{X) = 2.0 AND N= 1000 
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Table 3. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
WEIBULL MODEL WITH C, = 5.0, Q = 1.0, E{X,) = 2.0 AND N= 1000 
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Table 4. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
WEIBULL MODEL WITH C, = 8.0, Q = 1.0, £(^1) = 2.0 AND N= 1000 
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Table 6. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
WEIBULL MODEL WITH C, = 2.0, C, = 1.0, AND £(.\;) = 2.0 (N = 10) 
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Table 8. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
VVEIBULL MODEL WITH C, = 2.0. C, = 1.0, AND £(A;) = 2.0 (N = 250) 
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Figure 4. The Average Age Replacement Times For The Weibull Model 
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IV. SIMULATION SETTING FOR THE GAMMA DISTRIBUTION 


A. UNDERLYING LIFE DISTRIBUTION 

This chapter is concerned with the estimation of the optimal replacement time when 
it is known that the underlying lifetime distribution is a member of the two parameter 
Gamma family with shape parameter p > 1 and scale parameter 6 > 0, where the density 
is given by 

J{t) = — - for t>0. (4.1) 

•/w J 

lf/> = 1.0, the Gamma density with scale parameter 6 reduces to 

= i>0. (4.2) 

In equation (4.2), we have an exponential density with parameter 6. The reciprocal of 

the failure rate of the Gamma distribution is given in equation (4.3) 

-fi - \ffr' (4.3) 

By solving equation (4.3) analytically for different values of p, we get useful functional 
forms of r{i) in Table 9 on page 31. The failure rate for the Gamma random variable 
with ;»> 1.0 is a strictly increasing continuous function and is bounded above by 
A unique and finite optimum replacement policy cp* exists and will be finite if and only 
if (p — 1) is strictly greater than r ’ unplanned and planned 
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replacement costs (C, > Cj). Similar to the Weibull case, this unique minimum of C{(p) 
occurs at the minimum point v/heit the first derivative is zero. 


Table 9. FUNCTIONAL FORM OF rU) 


Shape param¬ 
eter p 

Scale parame¬ 
ter 6 

Functional form 

2.0 

1.0 

t 



t+ 1 

3.0 

0.666667 

27F 

18F+24r-i-16 

4.0 

0.5 

8/3 

4/3 + 6/3 -1- 6/ +3 

5.0 

0.4 

3125/" 

1310/" -}- 2000/3 + 2400/3 1920/ 4-768 


In this chapter cost C,, C, are chosen so that (p*is finite and unique. We used shape 
parameter p = 2.0. 3.0, 4.0, 5.0 (/’>!) and scale parameter 6 = 2 2, 2.3, 2/4, 2,5 in 
our simulation, respectively. This selection of p values gives us a range of distributions 
which become more like the exponential as p decreases from 5.0 to 2.0. To make fair 
comparisons between Gamma distributions, the scale parameter 9 was chosen so that 
the expected system lifetime L{X,) = 2.0. See Figures 6 and 7, for plots of the Gamma 
densities and corresponding failure rates. 
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GAMMA DENSITY FUNCTION 

f(t) = / r(p)9» 



Figure 6. The Gamma Density Function f(t) with E(Xi} = 2.0 
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GAMMA DENSITY FAILURE RATE FUNCTION 
lA(t) = 



Figure 7. The Failure Rate of The Gamma Distribution Mith £(A',) = 2.0 

B. OPTIMAL REPLACEMENT TIME 

When the distribution of a system lifetime is Gamma, the expected long run a\era'ge 
cost per unit time under a simple age replacement policy with scheduled replacement at 
age i. b> using equation (2.1) is given by 


C, 


n-\ 


riDO^ 


dx + Cj 


P-X 

x^ e 0 
V{p)0' 


dx 


Coji) = 


(4.4) 


/■(;r) dx 
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In equation (4.4), the corresponding survival function F{i) is, 


F{t) = 


n I i 

C 0 


r{p)(f 


dx, 


r>0. 


(4.5) 


In this thesis, we obtained the numerical values for F{t), by solving equation (4.5) ana¬ 
lytically for some specific values of p. Table 10 gives the useful functional form of 
F{i). See Figure 8 for plot of the Q,p(r) when the underlying life distribution is Gamma 
with shape parameter p varying from 2.0 to 5.0. For each curve on Figure 8 the optimal 
replacement time cp* can be located on the x-axis at the minimum point. 


Table 10. FUNCTIONAL FORM OF F(T) 


Shape param¬ 
eter p 

Scale parame¬ 
ter 0 

Functional form 

2.0 

1.0 

1+ 1 
e' 

3.0 

0.666667 

2.25F -r 3/ + 2 

2(?'-'' 

4.0 

0.5 

4F -f 6F 4- 6/ -f 3 

3c'»’- 

5.0 

0.4 

1310/^ -f- 2000r^ -1- 2400F -f 1920/ -f 768 

768(?=-^^ 
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LONG RUN EXPECTED AVERAGE COSTS FOR GAMMA DISTRIBUTION 


(Cj = 5.0 , C, = 1.0) 



Figure 8. The Long Run Expected Average Cost Cunes uith ^(A') = 2.0 


C. SEQUENTIAL ESTIMATION PROCEDURE OF THE GAMMA 
DISTRIBUTION. 

For the Gamma distribution, we will use our parametric sequential estimation pro¬ 
cedure again which we mentioned in Chapter 11, with the known shape parameter 
/> > 1.0 and the unknown scale parameter 0. 

/■ 

As in the Weibull case, minimizing Q (/) (where 0„ is the .MLE of d after n re- 
placements) can be simplified with the following result: 

Lemma 2. Let (p* minimize then cp* - where jt* minimizes C|_,(.v). 
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PROOF. From the cost function for the Gamma distribution in equation (4.4), 

JC 

after making the substitution y = ~^, the numerator may be written as 


’ er{p)d^ 


0 dy + 


f^OO - 1 « 

^ er{p)e^ 

J i 


6 dy. 


Canceling the scale parameters (0) we obtain 


(4.6) 


C, 


/»— 


/*oo 


^0 


rf/) 


dy + C2 


■>,-r 




dv. 


(4.7) 


The integrands of equation (4.7) are the density of the Gamma distribution with 6=1. 
Finally the numerator of (4.4) becomes 


C.r,=,(-^) + C2/-',=,(y). (4.8) 

In equation (4.4) the denominator may be written as 



dv du. 


(4.9) 


If we first solve the inner integral b> changing variable y with we get 
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f Av)(^v = 


OO V 


r{p)d^ 


p-\ -y 


r{p)e^ 


After cancellation, the denominator is 


n-1 _v 

I*" e ' 


(4.10) 


Again the integrand of equation (4.10) is the density of the Gamma distribution with 
6 = 1.. From (4.8) and (4.10), the cost function may be written as 


CeJj) = 


+ C,FUj) 


ro=iij)du 


(4.11) 


If we let the y = —, equation (4.11) becomes, 
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jc” 

Therefore, if;c* minimizes, and <p* minimizes Q_,(r) then (p* = —q~- ■ 

There is no closed form solution for the MLE of 0 based on right censored data. 
Thus, we use the EM algorithm to find the MLE numerically. 

1. EiM Algorithm 

" The Expeciaiion-Xfaxiinizaiion (EM) algorithm is an iterative algorithm used 
to compute the MLE in incomplete-data problems. The algorithm is applicable to more 
general missing-data patterns, but usually involves more computations than methods 
designed for special patterns. The iterations of the EM algorithm are 

1. Replace missing values by estimated values. 

2. Estimate parameters. 

3. Recstimate the missing values assuming the new parameter estimates are correct. 

4. Recstimate parameters. 

and so forth, iterating until convergence " (Ref 12: pp. 127-141). Since each iteration 
of the algorithm consist of an expectation step (E) followed by a maximization step (M) 
it calls the E.M algorithm. 

The M step is performed by maximizing the likelihood as if there were no miss¬ 
ing data. Thus, the .M step of E.M uses the identical computational methods as .M LE 
from/(0,pLV) (Ref. 12: pp. 127-141). With the shape parameter/> assumed known, the 
maximum likelihood estimator of 0 based on n uncensored observations X. ... X„ is 





( 4 . 13 ) 


The E step finds the conditional expectation of the missing (censored) data given 
the observed data and current estimated parameters, and then substitutes these expec¬ 
tations for the missing (censored) data [Ref. 12: pp. 127-141]. Thus at the Jth iteration 

A 

of the EM algorithm the .MLE 9„ based on n replacements is approximated by 


A 

= 


fi 


(4.14) 


where E{X\X > Z) is the conditional expectation of the random variable X ~ Ganmia( 
6 = ,_„/))• The functional forms of this conditional expectation for the specific shape 

parameter p are given in Table 11 on page 40. 

As in the Weibull case, we have no information about the age replacement time 
(p at the first obser\ation. Thus, the Z, value equals the first observation X^. We find 
(by using equation (4.13)). we get 


0, = IT- 


From the Lemma 2, we can estimate cp* such that where minimizes 


J^cp), and .V* minimizes Q ,,(.v). 

At the second obscrsalion. we have two cases. If Aj is less than (p*, then Oj is 
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Table 11. FUNCTIONAL FORM OF £(A:|X></>) WHERE X^GAMjO.p) 


Shape pa¬ 
rameter p 


Scale pa¬ 
rameter 6 


2.0 


1.0 


3.0 


0.666667 


4.0 


0.5 


5.0 


0.4 


(p' 


Functional form 


<jo^ -f- 2d(() -1- 26^ 

(p + 6 


(p^ + 3di(p^ + 2d(p + 26^) 
(p^ -1- 2B(p -1- 20^ 


(p^ + 40((p^ 30((p^ -H 2e<p + 26^)) 

(P^ + I0{(p'- + 20(p -f 20‘) 


-!- 59{(p^ 4ei(p^ -t- 30{(p^ + 2ecp -f- 2^-))) 
-f 40^^^ -1- 36(</>^ -f- 2d(p -}- 20')) 


O-i 


T| -f At 


Otherwise, the observation is censored, and we use the equation (4.14; for the E (ex¬ 
pectation) step of the EM algorithm until O 2 convergence. These iterations arc 


^2.1 ~ 


^2.2 ~ 


z, 4- i:;^L\\\x,>^p*-] 

2p 


Z, + U2lA2>^>n 





















and so on. Finally, these iterations converge 62 ^. When the difference of the absolute 

A A 

value of 02^, and 02^.i is small, the stopping criteria is satisfied. At that point, we can 

A A A A 

replace the 62 ^ with 62 - Again, we can estimate ^2 such that = - 7 —. 

A ^2 

The procedure is then repeated. After determining 6 values for each replace¬ 
ment, we can apply this estimated 6 values to the Lemma 2. The estimated optimal age 
replacement time can be expressed by. 


(4.15) 

en 

For large n, the estimated optimal age replacement time converges to an optimal age 
replacement time cp*. 

D. SIMULATION RESULTS FOR THE GAMIMA DISTRIBUTION 

1. Finding of The Scale Parameter 0 

The four different Gamma distributions used in the simulation have p values 2.0, 
3.0, 4.0. 5.0. and the scale parameter 6 was chosen so that the expected system lifetimes 
£(A’)= 2.0. The mean of the Gamma distribution is L(X) = p0. The 0 value can be ob¬ 
tain since both E(X) and p are known. For example, if /?= 3.0, then 0 = 0.666667. 

2. Finding of The Actual Optimal Age Replacement Time and Corresponding Actual 
Cost. 

Similar to the Weibull case, we found the minimum value of the cost function 
(4.4) and corresponding actual age replacement time (p* for different costs C,, Cj and 
parameters (0, p) by simulation. 

3. Finding of The Average Age Replacement Times and The Average Costs 

As in the Weibull case, in order to find the average of the estimated optimal age 
replacement times and the average actual costs based on lifetimes with a Gantma dis- 
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tribution, we wrote the Fortran simulation program Avergam in Appendix D. From 
Lemma 2, x* = (p*6 where the actual age replacement time is taken from the simu¬ 
lation program Sim and the 6 values are calculated from the equation E{X)=pd for the 
specific value of the shape parameter p. The unplanned and planned replacement costs 
C, and Cj (to calculate the cost function); the scale parameter 6 (to calculate the x* 
value); the actual age replacement time (from program Sim to calculate the x* value and 
to find MSE for each estimated age replacement time); the actual cost value (from pro¬ 
gram Sim to find MSE value for each estimated cost); and the shape parameter p must 
be given by the user in the initialization part of the simulation program. In much of the 
simulation, we considered the sample sizes A’ = 10, 50, 250 and 1000. The other sample 
sizes are also considered, but in much less detail. Each simulation is based on generating 
1000 sequences of A' sjstem lifetimes. We used these system lifetimes (A') one at a time 
for the sequential estimation procedure. At the first observation Z, = ^Vi, and d,= 1. 
As mentioned before, the scale parameter 0 values are estimated by the E.M algorithm. 
After finding 6, the estimated age replacement time values are calculated by using the 
equation (4.15) in the simulation. Similar to the Weibull case, by using A' generated 
system lifetimes, we determine A' estimated age replacement times and A' estimated costs 
in each simulation. We repeat this simulation 1000 times (NREP= 1000) and then we 
find the a\cragc \aluc for both age replacement times and replacement costs. See the 
equations (3.15) and (3.16) of Chapter 111. At each simulation we also calculated the 
MSE values for both the age replacement times and the long run expected average costs 
by using the equations (3.17) and (3.18). 

Tables 12, 13, 14 and 15 summarize the simulation results of the optimal age 
replacement times and the minimum long run expected replacement costs per unit time 
for different \aiues of p with A'= lOoO for fixed costs C, = 2.0, C, = 5.0, C, = 8.0. 
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C, = 10.0 while holding the Q fixed at 1.0. Tables 16, 17 and 18 summarize the simu¬ 
lation results using the three sample sizes of 10, 50 and 250 with different values of the 
shape parameter p for fixed costs C, = 2.0 and Q = 1.0. Included in Tables 12~18, is 
the probability that a system will fail before the time cp under the optimal age replace¬ 
ment policy (3.19). 

We have also plotted the results of the average age replacement times and the 
average costs for different shape parameter p. See Figure 9, for plot of the average age 
replacement times for p=3.0, 4.0, 5.0 {p=2.0 was not selected because its average age 
replacement times were high acording to the others) and Figure 10, for the average costs 
for p= 2.0, 3.0, 4.0, 5.0 when the fixed costs C^-2.0 and €2 = 1.0. 

From these plots and the results from Tables 12~18, we observe that when p 
decreases from 5.0 to 2.0 (i.e., the underlying life distribution is becomes more expo¬ 
nential) the long run expected optimal replacement costs increase. When the ratio of the 
unplanned replacement cost C, to the planned replacement cost Cj increases the optimal 
replacement time for a specific p decreases. The larger values of <p*, insure that a small 
percentage of replacement will be made before failure, which is what we desire if the 
system's life distribution is close to exponential. 

As in the Weibull case, by looking the tables for different sample sizes, we ob¬ 
serve that the average cost per unit time up to the A'th replacement decreases with .V, 
the number of replacement. Thi result is promising because the ultimate goal of the 
sequential estimation procedure is to decrease costs while sampling. Even though as 
A’ -+ 00 , the long run average cost per unit time will approach the optimal replacement 
cost C{(p*) with probability 1.0 [Ref. 11], there is no guarantee that the average re¬ 
placement cost will decrease for the first observations. For large sample sizes, our esti¬ 
mated average (p* and the aserage C(c*) values arc too close to their actual values. 
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Table 12. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
GAMMA MODEL WITH Q = 2.0, Q = 1.0, AND HA,) = 2.0 (N= 1000) 


Shape parameter p 

p=2.0 

II 

p = 4.0 

p=5.0 

Scale parameter 6 

I.O 

0.666667 

0.5 

0.4 

x* = 0x(p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal replacement time 
(P* 

35.35100 

3.17410 

2.20550 

1.89690 

Average 

35.37980 

3.17659 

2.20480 

1.89694 

MSEof<p* 

0.62197 

0.00371 

0.00168 

0.00106 

Long run expected optimal 
replacement cost C{(p*) 

1.00000 

0.99456 

0.97167 

0 93374 

Average 

l.OOOSl 

0.99552 

0.97201 

0.94614 

MSE of 

0.0(1049 

0.00036 

0.00032 

0.00041 

P( a; < (p*) 

1.00000 

0.61467 

0.18173 

0.03994 
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Table 13. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
GAMMA MODEL WITH C, = 5.0. C = 1.0. AND EiA,) = 2.0 (N= 1000) 


Shape parameter p 


II 

CO 

O 

p = 4.0 

p= 5.0 

Scale parameter 0 

1.0 

0.666667 

0.5 

0.4 

x* = 0x(p* 

1.30500 

0.67220 

0.48920 

0.39952 

Optimal replacement time 

<p* 

1.30500 

1.00830 

0.97840 

0.99880 

Average (p* 

1.30421 

1.00727 

0.97494 

0.99402 

MSE of (/>* 

0.00177 

0.00107 

0.00106 

0.00202 

Long run expected optimal 
replacement cost C{(p*) 

2.26480 

1.87695 

1.63237 

1.44237 

Average C(<p*) 

2.26830 

1.88581 

1.64236 

1.48797 

MSE of qc>*) 

0.00528 

0.00367 

0.O0307 

0.01029 

P>-V,<v>*) 

0.37495 

0.19428 

0.13517 

0.00291 


























































Table 14. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
GAMMA MODEL WITH C, = 8.0, Q = 1.0, AND ^A;.) = 2.0 (N= 1000) 


Shape parameter p 

p=2.0 

p=3.0 

II 

o 

p= 5.0 

Scale parameter B 

1.0 

0.666667 

0.5 

0.4 

x* = 0x<p* 

0.81890 

0.49540 

0.38620 

0.32956 

Optimal replacement time 

0.81890 

0.74310 

0.77240 

0.82390 

A^erage (p* 

0.81677 

0.73945 

0.76560 

0.81510 

MSEofcp* 

0.00124 

0.00112 

0.00172 

0.00334 

Long run expected optimal 
replacement cost 

3.15166 

2.28405 

1.97650 

1.68813 

Average C{<p*) 

3.16146 

2.10281 

2.00406 

1.77106 

MSE of C(<^*) 

0.01767 

0.01076 

0.01375 


P(A'. <<?>*) 

■i 

0.03960 

0.00SO6 

0.00120 


4(' 
























































Table 15. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 


GAMMA MODEL WITH C, = lO.O, C, = 1.0, AND £(A;) = 2.0 
(N= 1000) 


Shape parameter p 

p=2.0 

p=3.0 

tt 

o 

p = 5.0 

Scale parameter 6 

1.0 

0.666667 

0.5 

0.4 

X* = 6 X (p* 

0.68010 

0.43693 

0.35020 

0.30424 

Optimal replacement time 
<P* 

0.68010 

0.65540 

0.70040 

0.76060 

Average <p* 

0.67754 

0.65103 

0.68916 

0.74911 

MSE of <p* 

0.00113 

0.00111 

j 

0.00263 

0.00359 

Long run e.xpected optimal 
replacement cost 

3.64327 

2.64540 

2.14725 

1.80642 

Average C{<p*) 

3.66032 

2.67401 

2.19022 

1.90982 

MSE of C{<p*) 

0.02993 

0.01724 

0.02745 

0.04275 


P(.V,<(P*J 


892 


0.02894 


0.00576 


0.00082 
























































Table 16. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 


GAMMA MODEL WITH C, = 2.0, Q = 1.0, AND ^(A',) = 2.0 (N= 10) 


Shape parameter p 

p=2.0 

p=3.0 

P 

II 

cu 

II 

o 

Scale parameter B 

1.0 

0.666667 

0.5 

0.4 

x* = 0x(p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal replacement time 
(P* 

35.35100 

3.17410 

2.20550 

1.89690 

Average (p* 

37.57612 

3.23432 

2.18015 

1.82749 

MSE 

93.99036 

0.51523 

0.24643 

0.16467 

Long run expected optimal 
replacement cost C((p*) 

1.00000 

0.99456 

0.97167 

0.93374 

Average C{(p*) 

1.06294 

1.02846 

1.00494 

0.97408 

MSE of C{ip*) 

0.07521 

0.04433 

0.03426 

0.02609 


1.00000 

0.61467 

0.18173 

0.0399J 


4S 

























































Table 17. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
GAMMA MODEL WITH C, = 2.0. Q = 1.0. AND £(A;) = 2.0 (N = 50) 


Shape parameter p 

p=2.0 


II 

o 

p=5.0 i 

Scale parameter 6 

1.0 

0.666661 

0.5 

0.4 

1 

AT* = 0 X <p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal replacement time 

35.35100 

3.17410 

2.20550 

1.89690 

Average (p* 

35.84104 

3.19201 

2.20190 

1.87172 

MSE of <p* 

13.43787 

0.07847 

0.03622 

0.03245 

Long run expected optimal 
replacement cost C{(p*) 

1.00000 

0.99456 

0.97167 

0.93374 

Average C(<p*) 

1.01386 

! 

1.00393 

0.98294 

0.95588 

i 

MSE of C(<p*) 

0.01075 

0.00731 

0.00583 

0.00545 

P(X <<,'>*) 

1 

1.00000 

0.61467 

0.18173 

0.03094 
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Table 18. ESTIMATED OPTIMAL REPLACEMENT TIMES OF THE 
_GAMMA MODEL WITH C, = 2.0, C, = 1.0, AND £(A;) = 2.0 (N = 250) 


Shape parameter p 

P=2.0 

p=3.0 

q 

II 

a 

p = 5.0 

Scale parameter 0 

1.0 

0.666667 

0.5 

0.4 

x* = ex(p* 

35.35100 

2.11607 

1.10275 

0.75876 

Optimal replacement time 

(p* 

35.35100 

3.17410 

2.20550 

1.89690 

Average (p* 

35.47148 

3.18530 

2.20761 

1.89861 

MSE of (p* 

2.53150 

0.01452 

0.00649 

0.00436 

Long run expected optimal 
replacement cost C((p*) 

1.00000 

0.99456 

0.97167 

0.93374 

Average C{<p*) 

1.00340 

0.99883 

0.97499 

0.94996 

MSE of C((p*) 

0.00202 

0.00141 

0.00124 

0.00127 


1.00000 

0.61467 

0.18173 

0.03994 
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COUPABISON OF AGE REPLACEUENT TUIE 



Figure 9. The Average Age Replacement Times For The Gamma Model 


COMPARISON OF AVERAGE COSTS 
C,=2.0 , C,=1.0 



Figure 10. The Average Costs For The Gamma Model 
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V. CONCLUSIONS AND RECOMMENDATIONS 


Throughout this thesis, we have considered age replacement policies. In such poli¬ 
cies a unit is always replaced at the time of failure or units of time after its installation, 
whichever comes first. The time cp is called the scheduled replacement time or the 
scheduled censoring time. The optimal replacement time ip* achieves the minimum long 
run expected maintenance cost per unit time. The results of the simulations show that 
the parametric estimators work well under the conditions for which they were intended 
and by using the sequential estimation procedure significant cost and time savings can 
be cfiectcd. 

In most cases, the estimated optimal age replacement times and the actual costs got 
close to the true optimal age replacement time and the minimum expected cost per unit 
time respectively, even with moderate sample sizes. Therefore, we conclude that our 
sequential estimation procedure of the age replacement policies to minimize the long run 
expected cost can be applied to the real problem. 

An important part of our analysis, which would be very difficult to prove theore¬ 
tically. showed that the average actual cost per unit time decreased monotonically as the 
sample si/c increased. Such a decrease makes intuitive sense, i.c. an age replacement 
policy using an estimated (p* should do better as more data is collected. This result is 
desirable for a sequential estimation procedure. 

Directions for Future Research 

It is our hope the ideas and the sequential estimation procedure described in this 
thesis represent a plateau for the development of the optimum age replacement policies 
to minimi/e the long run expected maintenance costs. The minimization problem is in 




principle dilTicult and future developments will exploit fundamentally new ideas. Here, 
we briefly describe directions in which future research might be pursued. 

• We have described the scenario in which the underlying lifetimes are i.i.d., along 
with a straight forward age replacement policy. In many situations this model is 
hot adequate. For example, if an item is looaired at failure rather than replaced, 
the Lid. assumption is equivalent to requiring that the repaired item function as 
well as a new one. Clearly, modeling times between failure of a repairable system 
as Lid. is inappropriate. More realistic models incorporate the possibility that re¬ 
pairs are less than perfect. It is also possible that the quality of planned mainte¬ 
nance varies from time to time. If an imperfect repair model is permitted how 
should the sequential estimation procedure be changed to fit new situation. 

• Sequential estimation when the underlying life distribution F comes from the 
normal or the modified extreme value distributions which have increasing failure 
rate. 

• The minimization of the long run expected costs is not appropriate under all cir¬ 
cumstances because in this model the planned and unplanned costs arc fixed. 
Other cost functions need to be considered. For example, costs can be modeled to 
change with time. 

• How to change our estimation procedure if we have difl'erent systems such that 
serial, parallel or bridge systems with same or different lifetimes rather than one 
unit and one lifetime. 
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APPENDIX A. APL CODE TO CALCULATE LAMBDA OR MEAN FOR 

THE WEIBULL DISTRIBUTION 


V WEIBULLiAiEiLiX 

C13 fl THIS APL PROGRAM CALCULATES THE SCALE PARAMETER 
C23 fl LAMBDA AND THE EZXl FOR THE WEIBULL DISTRIBUTION. 
C33 R THE SHAPE PARAMETER a AND THE OTHER UNKNOWN 
[43 fl MUST BE GIVEN BY THE USER. 

[53 LUCK^O 

[63 START:'PLEASE TYPE 1 IF YOU WISH TO FIND LAMBDA' 

[73 'PLEASE TYPE 2 IF YOU WISH TO FIND EZXl' 

[83 TYPElf-Zi 

[93 -i-aTYPEl^'l' )AifYPEl*'2' ))/WARNINGl 

[103 ■>(1+2’YPS1= « 2 ' )/EXPECVA 

[113 LAMBDA:'PLEASE ENTER THE a VALUE' 

[123 M 

[133 'PLEASE ENTER THE EXPECTED VALUE EZXl' 

[143 M 

[153 fl LAMBDA = GAMMAil/a) / (ax Em) 

[163 X^UA 

[173 (I(X-l))+(vlxE) 

[183 ^END 

[193 EXPECVA:'PLEASE ENTER THE a VALUE' 

[203 

[213 'PLEASE ENTER THE LAMBDA VALUE' 

[223 M 

[233 fl Em = GAMMA il/a) / (ax LAMBDA) 

[243 X-e-l+i^ 

[253 (I (X-l))+(v?xL) 

[263 ^END 
[273 WARNINGl: 
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[28] LUCK<-LUCK+1 

[29] -^iLUCK^D/CHEWl 

[30] 'YOU CRN ONLY ENTER THE NUMBERS 1 OR 2' 

[31] ^START 

[32] CHEWli'JUST TYPE 1 OR 2' 

[33] ^START 

[34] WARNING!', 

[3 5] LUCK<-LUCK+1 

[36] ■i-UUCK^2)/CHEW2 

[37] 'YOU CAN ONLY ENTER THE LETTERS Y OR N' 

[38] ^END 

[39] CHEW2:'JUST TYPE Y OR N' 

[UO] -^END 

[41] END:'DO YOU WANT TO RUN THE PROGRAM AGAIN?(Y/N)' 
[4 2] TYPE2^\Ii 

[43] ^(il^TYPE2)='Y')/START 

[44] ^(.l^TYPEl^'Y' ')f^il^TYPE2xN)/WARNING2 

[45] 'BYE BYE' 

V 






APPENDIX B. APL CC TO CALCULATE THE ACTUAL AGE 
REPLACEMENT TIME AND CORRESPONDING THE ACTUAL 

MINIMUM COST 


V SmiCliC2iIiJiFXiXiT;XAiYAiCiDiXMINiYMIN 
CID P THIS PROGRAM SIMULATES THE COST FUNCTION (EQ. 2.1) TO 

C2] R FIND THE MINIMUM VALUE (YMIN) OF THE COST FUNCTION AND 

C3] R CORRESPONDING AGE REPLACEMENT TIME (XMIN) FOR THAT 

[4] R POINT. AFTER FINDING MINIMUM VALUES INSIDE THE LOOPl 
[53 R ir REPEATS THE PROCEDURE 300 TIMES INSIDE THE L00P2. 
[63 R FINALLY, THE PROGRAM GIVES THE AVERAGE VALUES FOR 

[73 R BOTH MINIMUM POINT AS AXST AND ACST. 

[83 R 

[93 ff-f-CiSOOO + lOO 

[lO: R THIS GIVES US A VECTOR OF T{0.01, 0.02, .... 50) 

[113 R TO CALCULATE FIRST C(O.Ol) AND THEN £7(0.02) UP TO 
[123 R £7(50) OF 5000 COST VECTOR. 

[133 R INITIALIZATION... 

[143 fi UNPLANNED AND PLANNED REPLACEMENT COSTS MUST BE 
[153 fl GIVEN BY THE USER. 

[163 £71-t-5 

[173 £72-«-l 

[183 M^iO 
[193 W-f-iO 
[203 J^O 

[213 R J IS THE INCREMENT OF THE LOOP2 J-1, 2, ..., 300 
[223 R MODEL... 

[233 L00P2: 

[243 X-(-50G0 WEIRAND 2 2.2567587 

[253 fl PREVIOUS LINE, GENERATES 5000 SYSTEM LIFETIMES FROM 
[263 fl WEI{ALPHA=2.0,BETA=2.2501331) AS VECTOR X.HERE BETA 
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[27] p VALUE REPRESENTS 1 OVER LAMBDA=1*0.44311346 

[28] n FOR GAMMA DISTRIBUTION LINE 24 CAN BE SWITCH WITH 

[29] R Xf-5000 GAMRAND 4 0.5 GAM( P=4 tTHETA-0,5 ). 

[30] JfrJ+1 

[31] C-e-iO 

[32] I-c-0 

[33] n I IS THE INCREMENT OF THE INNER LOOP 1=1, 2, 5000 

[34] LOOPli 

[35] I^I+l 

[36] R C IS THE SIMULATED COST FUNCTION 

[37] R D IS THE DENUMERATOR OF THE COST FUNCTION 

[38] D-f-((+/(XL2’[J])) + 5000) 

[39] Cf-C, ( ( (C2X (1-FX) ) + (Clx (FX-«-( (+/X5Z’[I] )+5000 ) ) ) )*D) 

[40] R IN THE FIRST LOOP C VECTORS OBTAIN FOR EACH T 

[41] •>(I<5000)/LOOP1 

[42] YMIN-^l/C 

[43] XMINt-TZli-iCl 

[44] R YMIN: THE MINIMUM VALr’E OF THE COST FUNCTION 

[45] R FOR SPECIFIC T 

[46] R XMINi THE CORRESPONDING AGE REPLACEMENTTIME (I) 

[47] XA<rXA,XMIN 

[48] YA<-YA,YMIN 

[49] R XAi THE VECTOR OF THE AGE REPLACEMENT TIMES (300) 

[50] fl YAi THE VECTOR OF THE YMIN (300) 

[51] ■*-iJ<20Q)/L00P2 

[52] AXST^i+/XA)*QXA 

[53] ACST*-i+/YA')*pYA 

[54] R AXST: THE AVERAGE VALUE OF THE AGE REPLACEMENT 

[55] R TIMES AFTER 300 REPETITIONS. 

[56] fl ACST: THE AVERAGE VALUE OF THE YMIN 

[57] fl AFTER 300 REPETITIONS. 

V 
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APPENDIX C PROGRAM AVERWEIB 


c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM AVERWEIB 


Wf* ****** 


*** PURPOSE ; This program calculates the average costs and *** 
*** corresponding average replacement times by using the *** 
*^c* sequential estimating procedure for the weibull distri- *** 
*** bution. The program also calculates the mean square *** 
*** error values for the average costs and the average rep- *** 
*-f* lacement times at the each run. *** 


PARAMETER (N=1000, NREP=1000) 

INTEGER I, J, K, L, M, SEED, DEL(N) 

REAL*8 X(N), LAMHAT(N), Z(N), FISTAR(N), XSTAR, SUMZ(N), NUM(N) 
REAL*8 SUM(N), SUM2(N), AVG(N), FSTAR, MSEAGE(N), Cl, C2, COSTST 
REAL*8 DEN(N), ACOST(N), SUM3(N), SUM4(N), AVACOS(N), MSECOS(N) 


c 

*** 

The following real 

numbers are not double precision 

ieicit 

c 


because the random 

number generator subroutine usee 

*** 

c 

*** 

single precision 

only. 

*** 


REAL U,LAM,AL 




c 

***y 

.'**************iV*- 

Tin 

* Vc Vir #V Vc Vc 

TiciciTic 

c 

Vr** 

KEY VARIABLES 



icici: 

c 


ISEED 


The random number seed 

*** 

c 


X 


Generated system lifetimes 

**yr 

c 

*** 

LAM 


The scale parameter 

Vf^ryc 

c 

*** 

AL 


The shape parameter 

iricic 

c 


Cl 


Unplanned replacement cost 

iciriT 

c 


C2 


Planned replacement cost 

icicif 

c 

*** 

XSTAR 


A constant which minimizes the cost 

iritic 

c 

*** 



function when the scale parameter 

*y«v 

c 

*** 



equals 1 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FSTAR 


COSTS! 


ACOST 

FISTAR 


LAMHAT 


MSEAGE 


AVACOS 


MSECOS 


The actual age replacement time (from *** 
the APL program Sim) *** 
The actuai I'ost (from the APL program *** 
Sim) *** 
The return generated random number *** 
The number of the run *** 
The number of the repetition *** 
Indicates the summation of the delta *** 
values to calculate the average costs *** 
The numerator of the cost function *** 
which given in the Equation (2.6) *** 
The denominator of the cost function *** 
which given in the Equation (2.7) *** 
The average cost *** 
The average age replacement time *** 
The minimum value of the X(i) (genera-*** 
ted system lifetimes) or the age rep- *** 
lacement time *** 
The estimated scale parameter *** 
The average value of the age replace- *** 
ment times after NREP repetition *** 
Mean Square Error for the average of *** 
the age replacement times *** 
The average value of the costs after *** 
NREP repetition *** 
Mean Square Error for the average of *** 
the costs *** 


*** INITIALIZATION *** 

«V yu 

SEED=16807 
LAM=0. 4557111655 
AL=1. ADO 


Cl=10. ODO 
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c 


c 

c 

c 


C2=l. ODO 
FSTAR=0.91690D0 

*vov From the main lemma xstar=fstar/lambda 

XSTAR=FSTAR/LAM 

C0STST=4. 048919D0 


DO 2 1=1,N 

SUM(I)=0. ODO 
AVG(I)=0. ODO 
MSE(I)=0.ODO 
SUM2(I)=0. ODO 
SUM3(I)=0. ODO 
SUM4(I)=0. ODO 
CONTINUE 


***Vf********************************************************ff*** 


*** CALCULATIONS 






DO 3 L=1,NREP 
C 

CALL EXCMS ('FILEDEF 11 DISK ALU OUTPUT A') 

CALL RANNUM (1,SEED,0. 0,1. 0,0,U) 

C 

X(!)=((-1*L0G(1-U))/(LAM**AL))**(1/AL) 

Z(1)=X(1) 

LAMHAT(1)=1/X(1) 

FISTAR(1)=LAMHAT(1)*XSTAR 
SUMZ(1)=X(1)**AL 

C *** First delta value is 1 (Failure occur before the time *** 

C *** and at the first observation there is no comparison for *** 

C *** the miniraums. Z(l) = X (1) 

DEL(1)=1 

C *** Calculation of the average cost for the first observation*** 

K’T’V'' 1 N—1 

noav, i j—ui 

DEN(1)=Z(1) 

AC0ST(1)=NUM(1)/DEN(1) 
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CALL RANNUM (1,SEED,0.0,1.0,0,U) 


X(J)=((-l*LOG(1-U))/(LAM**AL))**(1/AL) 

*** Comparisons to find the Z and delta values 
IF(FISTAR(J-1) .LT. X(J)) THEN 

SUMZ(J)=SUMZ(J-1)+(FISTAR(J-1)**AL) 
DEL(J)=DEL(J-1) 

Z(J)=FISTAR(J-1) 

ELSE 

SL'MZ(J)=SUMZ(J-1)+(X(J)**AL) 

DEL(J)=DEL(J-1)+1 

Z(J)=X(J) 

END IF 

LAMHAT(J)=(DEL(J)/SUMZ(J))**(1/AL) 

FISTAR(J)=LAMHAT(J)*XSTAR 

NUM(J)=(C1*DEL(J))+(C2*(J-DEL(J))) 

DENXJ)=DEN’(J-1)+Z(J) 

AC0ST( J) =N'UM( J ) /DEN( J ) 

CONTINUE 
DO 10 K=1,N 

SUM(K)=SUM(K)+FISTAR(K) 

SUM2(K)=SUM2(K)+(FISTAR(K)-FSTAR)**2 
SUM3(K)=SUM3(K)+ACOST(K) 
SUM4(K)=SUM4(K)+(AC0ST(K)-C0STST)**2 
CONTINUE 
CONTINUE 
DO 20 1=1,N 

AVG(I)=SUM(I)/NREP 
MSEAGECI)=SUM2(I)/NREP 
AVACOSCI)=SUM3(I)/NREP 
MSECOSCI)=SUM4(I)/NREP 

WRITE(11,*) AVG(I),AVACOS(I),MSEAGE(I),MSECOS(I) 






20 CONTINUE 
STOP 
END 

SUBROUTINE RANNUM(DISTN, SEED, RPARMl, RPARM2, IPARM, X) 

C *** this SUBROUTINE PROVIDES AN INTERFACE WITH THE LLRANDOMII*** 

C *** ROUTINES PROVIDED IN THE NONIMSL LIBRARY. THE PARAMETER *** 

C *** REQUIREMENTS AND CALLING PROCEDURES ARE AS FOLLOWS: *** 

Q 'jVfr* icic-k 

c *** dISTN = DISTRIBUTION TYPE YOU WANT TO SELECT *** 

C *** AN INTEGER BETW'EEN 1 AND 7 *** 

C *** SEED = THE RANDOM NUMBER SEED YOU WISH TO USE *** 

C *** RPARMl, RPARM2, AND IPARM ARE REAL AND INTEGER PARAMETERS*** 

C *** PASSED TO THE ROUTINE WITH MEANINGS WHICH VARY WITH THE *** 

C *** type OF DISTRIBUTION YOU DESIRE *** 

C *** X = THE RETURNED RANDONM NUMBER, IT IS ALWAYS REAL *** 

C *** 

C *** DISTRIBUTION NUMBERS AND THE ASSOCIATED FARM DEFINITIONS *** 

Q *** 

c *** 1--UNIFORM ON THE INTERVAL RPARMl TO RPARM2 *** 

C *** 2--NORMAL WITH MEAN RPARMl AND VARIANCE RPARM2 *** 

C *** 3--EXPONENTIAL WITH RATE RPARMl 

C *** 4--C0UCHY WITH A = RPARMl AND B = RPARM2 *** 

C *** 5--GAMMA WITH SrlAPE RPARM2 AND RATE RPARMl *** 

C *** 6--POISSON WITH RATE RPARMl *** 

C *** 7--GEOMETRIC WITH P = RPARMl *** 

REAL RPARMl, RPARM2, X 
I.NTEGER DISTN, SEED, IPARM, N 

Q **5V *** 

REAL TEMP, VARIAT(l) 

IF (DISTN. LE.O. OR. DISTN. GT. 8) THEN 

WRITEdO, *) 'ILLEGAL CALL TO RANNUM, BAD DISTN' 
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STOP 

ENDIF 

C *** 

GOTO (10, 20, 30, 40, 50, 60, 70), DISTN 
C *** *** 

c *** generate a uniform between RPARMl AND RPARM2 *** 

10 CONTINUE 

C *** *** 


IF (RPARMl - RPARM2.EQ. 0) THEN 

WRITEdO, *) 'ILLEGAL EQUAL RPARMS IN RANNUM' 

STOP 

ENDIF 

IF (RPARMl.GT.RPARM2) THEN 
TEMP = RPARMl 
RPARMl = RPARM2 
RPARM2 = TEMP 
ENDIF 

CALL LRND(SEED, VARIAT, 1, 1, 0) 

VARIAT(l) = RPARMl + (RPARM2 - RPARMl) * VARIAT(l) 

GOTO 99 

C GENERATE A NORMAL WITH MEAN RPARMl AND STDDEV RPARM2 *** 

20 CALL LNORM(SEED, VARIAT, 1, 1, 0) 

WRITE(*, *) 'NORMAL (0, 1) ', VARIAT(l) 

VARIAT(l) = (VARIAT(l) * RPARM2) + RPARMl 
GOTO 99 

Q *** 

C *** GENERATE AN EXPONENTIAL WITH RATE (1/MEAN) RPARMl *** 

30 CONTINUE 

IF (RPARMl.EQ.O) THEN 

WRITEdO, *) 'ILLEGAL ZERO RATE IN RANNUM' 

STOP 

ENDIF 

CALL LEXFN(SEED, VARIAT, 1, 1, 0) 


63 







VARIAT(l) = VARIAT(l) / RPARMl 
GOTO 99 

Q *** *** 

C ***GENERATE A COUCHY WITH A = RPARMl AND B = RPARM2 *** 

40 CONTINUE 

IF (RPARM2.LE, 0) THEN 

WRITEdO, *) 'ILLEGAL COUCHY SPREAD IN RANNUM, B = ' ,RPARM2 
STOP 
ENDIF 

CALL LCCHY(SEED, VARIAT, 1, 1, 0) 

VARIAT(l) = (VARIAT(l) * RPARM2) + RPARMl 
GOTO 99 

50 CONTINUE 

IF (RPARMl. LE.O) THEN 

WRITEdO, *) 'ILLEGAL NONPOSITIVE GAMMA RATE IN RANNUM' 

STOP ■ 

ENDIF 

IF (RPARM2.LE. 0) THEN 

WRITEdO, *) 'ILLEGAL SHAPE PARAMETER IN RANNUM' 

STOP 

ENDIF 

CALL LGAMACSEED, VARIAT, 1, 1, 0, RPARM2) 

VARIATd) = VARIATCl) * (1.0 / RPARMl) 

GOTO 99 

60 CONTINUE 

IF (RPARMl. LE.O) THEN 

WRITEdO, *) 'ILLEGAL POISSON RATE IN RANNUM' 

STOP 

ENDIF 

CALL LPOIS(SEED, VARIAT, 1, 1, 0, RPARMl) 

GOTO 99 

70 CONTINUE 

IF (RPARMl. LE.O) THEN 

WRITEdO, *) 'ILLEGAL GECM PROB IN RANNUM' 
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STOP 

ENDIF 

CALL LGEOM(SEED, VARIAT, 1, 1, 0, RPARMl) 

GOTO 99 

*** *** 
CONTINUE 
X = VARIAT(l) 

END 


65 



APPENDIX D. PROGRAM AVERGAM 


c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM AVERGAM 


*** PURPOSE : This program calculates the average costs and *** 
*** corresponding average replacement times by using the *** 
*** sequential estimating procedure for the gamma distribu- *** 
*** tion. The program also calculates the mean square error*** 
*** values for the average costs and the average replacement*** 
*** times at the each run. *** 

**5V****5V******?V***5V****Vr*************************A************** 


PARAMETER (N=1000, NREP=1000) 

INTEGER I, J, K, L, M, ISEED, DEL(N), DELTA(N) 

REAL*8 X(N), TEHAT(N), Z(N), FISTAR(N), XSTAR, COSTST 
REAL*8 SUM(N), SUM2(N), AVG(N), FSTAR, MSEAGE(N), Cl, C2 
REAL*8 ACOST(N), SUM3(N), SUM4(N), AVACOS(N), MSECOS(N) 
REAL*8 NUMCO(N), DENCO(N), Tl, T, NUM, A, B, C, D 


c 


The following Real 

numbers are not Double precision 


c 

*** the number generator subroutine uses single precision 

REAL U,TETA,P,RATE 

*** 

c 

■sV'jV* 

KEY VARIABLES : 


Vf7V9V 

c 

Vf Vf Vv 

ISEED ; 

The random number seed 

*** 

c 

*** 

X : 

Generated system lifetimes 

VfiV* 

c 


TETA : 

The scale parameter 


c 

*** 

P : 

The shape parameter 

*** 

c 


RATE : 

Reciprocal of the scale parameter to 


c 

*** 


use the subroutine Rannum 

•k'h'k 

c 

*** 

Cl : 

Unplanned replacement cost 

*** 

c 

*** 

C2 : 

Planned replacement cost 


c 

*** 

XSTAR : 

A constant which minimizes the cost 


c 



function when the scale parameter 

VfVf* 

c 

*** 


equals 1 

icici: 
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c 

*** 

FSTAR 

: The actual age replacement time (from 


c 

*** 


the APL program Sim) 

*** 

c 


COSTS! 

: The actual cost (from the APL program 


c 

*** 


Sim) 

*** 

c 

*** 

U 

: The return generated random number 


c 

*** 

N 

: The number of the run 


c 

*** 

NREP 

: The number of the repetition 


c 

*** 

DELTA 

: Indicates 0 or 1. If failure occurs 


c 

*** 


before the time age, then equals 0. 

*** 

c 

*** 


Otherwise equals 1 


c 

icicfc 

DEL 

: Indicates the summation of the delta 


c 

*** 


values to calculate the average costs 

•h-kic 

c 

*** 

NUMCO 

: The numerator of the cost function 


c 



which given in the Equation (2.6) 

ickk 

c 

Vr** 

DENCO 

: The denominator of the cost function 

*** 

c 

*** 


which given in the Equation (2.7) 

*** 

c 

*** 

ACOST 

: The average cost 


c 


FISTAR 

: The average age replacement time 

kkk 

c 

^VVcVc 

Z 

: The minimum value of the X(i) (genera 


c 



ted system lifetimes) or the age rep¬ 


c 

Vr*^V 


lacement time 

kkk 

c 

*yf* 

TEHAT 

; The estimated scale parameter 

kkk 

c 

Vf** 

T1 

: The converged scale parameter value 

kkk 

c 



at the E step of the EM algorithm 


c 


T 

: The value of the scale parameter at 

kkk 

c 

*** 


the previous calculation 

kkk 

c 


NUM 

: The numerator of the Equation 4. 14 

kkk 

c 


A, B, C, D 

: To determine the conditional expecta¬ 

kkk 

c 

vV** 


tion value from Table 8 

kkk 

c 

*** 

AVG 

; The average value of the age replace¬ 

kkk 

c 



ment times after NREP repetition 

kkk 

c 


MSEAGE 

: Mean Square Error for the average of 

kkk 

c 

>V>VVf 


the age replacement times 

kkk 

c 


AVACOS 

: The average value of the costs after 
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c 

c 

c 

c 

c 

c 


c 


MSECOS 


NREP repetition *** 

: Mean Square Error for the average of *** 
*** the costs 

*** TMTTTAT.T7ATTnM * 


ISEED=16807 
TETA=2.0/5. 0 
RATE=1. 0/TETA 
P=5.ODO 
Cl=2.ODO 
C2=l.ODO 
FSTAR=1. 89690D0 

*** From the main lemma xstar=fstar*theta 


XSTAR=FSTAR*TETA 
C0STST=0. 933742D0 
DO 2 1=1,N 

SUM(I)=0. ODO 
AVG(I)=0. ODO 
MSE(I)=0.ODO 
SUM2(I)=0. ODO 
SUM3(I)=0.0D0 
SUM4(I)=0. ODO 
2 CONTINUE 

Q Vf * * * * Vf “V i: Vc Vr Vr Vf Vr ■/: V: Vr Vf Vf Vr Vf * V? it itit it it itit it it itititic it it itititii it it ititititit itit it itit it it it ititit itititit 

C *** CALCULATIONS *** 

Q it it ititititit it it it itit it itit itit it itit itititit itit itit itit itit ititit it it it it itit itit it it it itit itit itititititititit itit itititit 


DO 3 M=1,NREP 
C 

CALL EXCMS ('FILEDEF 15 DISK P5C2 OUTPUT A’) 
CALL RANNUH (5,ISEED,RATE,P,0,U) 

C 

X(1)=DBLE(U) 
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Z(l)=x(l) 

C *** First theta value from the Maximum Likelihood Estimation *** 
TEHAT(1)=Z(1)/P 

C *** Use the main lemma to find the estimated age replac. time*** 

FISTARC1)=XSTAR/TEHAT(1) 

C *** First DELTA value is 1 (Failure occur before time age) *** 

C *** and at the first observation there is no comparison *** 

C *** Z(1)=X(1) *** 

DEL(1)=1 
DELTA(1)=1 

C *** Calculation of the average cost for the first observation*** 
NUMC0(1)=C1 
DENC0(1)=Z(1) 

ACOST(1)=NUMCO(1)/DENCO(1) 

DO 5 J=2,N 
C 

CALL RANNUM (5,1SEED,RATE,P,0,U) 

C 

X(J)=DBLE(U) 

C *** Comparison to find the Z and DELTA values *** 

IF(FISTAR(J-1) .LT. X(J)) THEN 
Z(J)=FISTAR(J-1) 

DEL(J)=DEL(J-1) 

DELTA(J)=0 

ELSE 

Z(J)=X(J) 

DEL(J)=DEL(J-1)+1 
DELTA(J)=1 
END IF 
T=TEHAT(J-1) 

90 CONTINUE 

NUM=Z(1) 

DO 100 L=2,J 
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IF(DELTA(L) .EQ. 1) THEN 
C *** Observation is uncensored 

NUM=NUM+Z(L) 

ELSE 

C *** Observation is censored 

C *** Calculate the conditional Expectation E(XlX>Z) values *** 
C *** by using Table 8 

A=(Z(L)*Z(L))+(2*Z(L)*T)+(2*T*T) 

B=( Z( L)*Z( L)*Z( L) )+( S^T’VA) 

C=(Z(L)*Z(L)*Z(L)*Z(L))+(4*T*B) 

D=( Z( L)*Z( L)*Z( L)*Z( L) )+( 5*T^^C ) 

NUM=NUM+(D/C) 

END IF 

100 CONTINUE 

T1 = NUM / (J*P) 

IF (ABS(Tl-T) .LT. 0.001) GO TO 200 

C *** This if statement is to check stopping criteria *** 

C *** If this satisfied we can accept T1 converged to the theta*** 

T=T1 

GO TO 90 
200 CONTINUE 

TEHAT(J)=T1 

FISTAR(J)=XSTAR/TEHAT(J) 

NUMC0(J)=(C1*DEL(J))+(C2*(J-DEL(J))) 

DENC0(J)=DENC0(J-1)+Z(J) 

AC0ST(J)=NUMC0(J)/DENC0(J) 

5 CONTINUE 

DO 10 K=1,N 

SUM(K)=SUM(K)+FISTAR(K) 
SUM2(K)=SUM2(K)+(FISTAR(K)-FSTAR)**2 
SUMS(K)=SUM3(K)+AC0ST(K) 
SUM4(K)=SUM4(K)+(AC0ST(K)-C0STST)**2 
10 CONTINUE 

3 CONTINUE 
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DO 30 1=1,N 

AVG(I)=SUM(I)/NREP 
MSEAGE(I)=SUM2(I)/NREP 
AVACOS(I)=SUM3(I)/NREP 
MSECOSCI)=SUM4(I)/NREP 

WRITE(15,300) AVG(I),AVACOS(I),MSEAGE(I),MSECOS(I) 
300 FORMAT (FIO. 7,' ',F10.7,' ',F10.7,' 

30 CONTINUE 
STOP 
END 




' ,F10. 7) 
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